ASV level
level="ASV"
path_maaslin="../intermediate_files/maaslin/Q1_norway/ASV/"
raw_linda_results[[segment]] <- list()
linda_results[[segment]] <- list()
supplements_wb <- createWorkbook()
pre_ltx vs healthy
group <- c("healthy","pre_ltx")
comparison_name <- paste0(group[1], " vs ",group[2])
linDA
# prepare the data
linda_data <- binomial_prep(colon_asv_tab,
colon_taxa_tab,
colon_metadata,
group, usage="linDA")
## Removing 900 ASV(s)
## Removing 78 ASV(s)
filt_colon_uni_data <- linda_data[[1]]
filt_colon_uni_taxa <- linda_data[[2]]
filt_colon_uni_metadata <- linda_data[[3]]
# fit the model
linda.obj <- linda(filt_colon_uni_data,
filt_colon_uni_metadata,
formula = '~ Group + (1|Patient)',prev.filter = 0.2)
## 339 features are filtered!
## The filtered data has 273 samples and 174 features will be tested!
## Pseudo-count approach is used.
## Fit linear mixed effects models ...
## Completed.
linda.output <- linda.obj$output
linda.output <- linda_renaming(linda.output, group)
group1 <- paste0(group[1], " vs ","Group",group[2])
for (grp in c(group1)){
raw_linda_results[[segment]][[grp]] <-
rawlinda.df(linda.output,
grp,
filt_colon_uni_data,
filt_colon_uni_taxa)
linda_results[[segment]][[grp]] <-
linda.df(linda.output,
grp,
filt_colon_uni_data,
filt_colon_uni_taxa)
}
# volcano plot
volcano_1 <- volcano_plot_linda(linda.output, group1,
taxa_table = filt_colon_uni_taxa) +
ggtitle(paste(group,collapse=" vs "))
volcano <- ggarrange(volcano_1)
# see the plot
volcano

MaAsLin2
Volcano plot
volcano1 <- volcano_plot_maaslin(fit_data,filt_colon_uni_taxa) +
ggtitle(paste(group[1], "vs", group[2]))
volcano <- ggarrange(volcano1)
# see the results
volcano
## Warning: ggrepel: 120 unlabeled data points (too many overlaps). Consider increasing max.overlaps

Group - Intersection
intersection_results <- group_intersection(group, list_intersections, list_venns,
linda.output, fit_data,
raw_linda_results,
segment = segment,
level=level)
list_intersections <- intersection_results[[1]]
list_venns <- intersection_results[[2]]
venn <- intersection_results[[3]]
# show the results
venn

Basic statistics
uni_df <- merge(basic_univariate_statistics(linda_data,group),
raw_linda_results[[segment]][[group1]],
by="SeqID",all=TRUE)
## [1] "healthy"
## [1] "pre_ltx"
## [1] "healthy"
## [1] "pre_ltx"
uni_df[["final_sig"]] <- uni_df$SeqID %in% list_intersections[[paste(segment,level,comparison_name)]][["SeqID"]]
uni_statistics[[segment]][[paste(level,comparison_name)]] <- uni_df
# for comparison
new_name <- paste(level,comparison_name)
addWorksheet(supplements_wb, sheetName = new_name)
writeData(supplements_wb, sheet = new_name, uni_df, rowNames=FALSE)
pre_ltx vs post_ltx
group <- c("pre_ltx","post_ltx")
comparison_name <- paste0(group[1], " vs ",group[2])
linDA
# prepare the data
linda_data <- binomial_prep(colon_asv_tab,
colon_taxa_tab,
colon_metadata,
group, usage="linDA")
## Removing 848 ASV(s)
## Removing 45 ASV(s)
filt_colon_uni_data <- linda_data[[1]]
filt_colon_uni_taxa <- linda_data[[2]]
filt_colon_uni_metadata <- linda_data[[3]]
# fit the model
linda.obj <- linda(filt_colon_uni_data,
filt_colon_uni_metadata,
formula = '~ Group + (1|Patient)')
## 0 features are filtered!
## The filtered data has 349 samples and 304 features will be tested!
## Pseudo-count approach is used.
## Fit linear mixed effects models ...
## Completed.
linda.output <- linda.obj$output
linda.output <- linda_renaming(linda.output, group)
group1 <- paste0(group[1], " vs ","Group",group[2])
for (grp in c(group1)){
raw_linda_results[[segment]][[grp]] <-
rawlinda.df(linda.output,
grp,
filt_colon_uni_data,
filt_colon_uni_taxa)
linda_results[[segment]][[grp]] <-
linda.df(linda.output,
grp,
filt_colon_uni_data,
filt_colon_uni_taxa)
}
# volcano plot
volcano_1 <- volcano_plot_linda(linda.output, group1,
taxa_table = filt_colon_uni_taxa) +
ggtitle(paste(group,collapse=" vs "))
volcano <- ggarrange(volcano_1)
# see the plot
volcano

MaAsLin2
volcano1 <- volcano_plot_maaslin(fit_data,filt_colon_uni_taxa) +
ggtitle(paste(group[1], "vs", group[2]))
volcano <- ggarrange(volcano1)
# see the results
volcano

Group - Intersection
intersection_results <- group_intersection(group, list_intersections, list_venns,
linda.output, fit_data,
raw_linda_results,
segment=segment,
level=level)
list_intersections <- intersection_results[[1]]
list_venns <- intersection_results[[2]]
venn <- intersection_results[[3]]
# show the results
venn

Basic statistics
uni_df <- merge(basic_univariate_statistics(linda_data,group),
raw_linda_results[[segment]][[group1]],
by="SeqID",all=TRUE)
## [1] "pre_ltx"
## [1] "post_ltx"
## [1] "pre_ltx"
## [1] "post_ltx"
uni_df[["final_sig"]] <- uni_df$SeqID %in% list_intersections[[paste(segment,level,comparison_name)]][["SeqID"]]
uni_statistics[[segment]][[paste(level,comparison_name)]] <- uni_df
# for comparison
new_name <- paste(level,comparison_name)
addWorksheet(supplements_wb, sheetName = new_name)
writeData(supplements_wb, sheet = new_name, uni_df, rowNames=FALSE)
post_ltx vs healthy
group <- c("healthy","post_ltx")
comparison_name <- paste0(group[1], " vs ",group[2])
linDA
# prepare the data
linda_data <- binomial_prep(colon_asv_tab,
colon_taxa_tab,
colon_metadata,
group, usage="linDA")
## Removing 1180 ASV(s)
## Removing 72 ASV(s)
filt_colon_uni_data <- linda_data[[1]]
filt_colon_uni_taxa <- linda_data[[2]]
filt_colon_uni_metadata <- linda_data[[3]]
# fit the model
linda.obj <- linda(filt_colon_uni_data,
filt_colon_uni_metadata,
formula = '~ Group + (1|Patient)')
## 0 features are filtered!
## The filtered data has 207 samples and 528 features will be tested!
## Imputation approach is used.
## Fit linear mixed effects models ...
## Completed.
linda.output <- linda.obj$output
linda.output <- linda_renaming(linda.output, group)
group1 <- paste0(group[1], " vs ","Group",group[2])
for (grp in c(group1)){
raw_linda_results[[segment]][[grp]] <-
rawlinda.df(linda.output,
grp,
filt_colon_uni_data,
filt_colon_uni_taxa)
linda_results[[segment]][[grp]] <-
linda.df(linda.output,
grp,
filt_colon_uni_data,
filt_colon_uni_taxa)
}
# volcano plot
volcano_1 <- volcano_plot_linda(linda.output, group1,
taxa_table = filt_colon_uni_taxa) +
ggtitle(paste(group,collapse=" vs "))
volcano <- ggarrange(volcano_1)
# see the plot
volcano
## Warning: ggrepel: 14 unlabeled data points (too many overlaps). Consider increasing max.overlaps

MaAsLin2
volcano1 <- volcano_plot_maaslin(fit_data,filt_colon_uni_taxa) +
ggtitle(paste(group[1], "vs", group[2]))
volcano <- ggarrange(volcano1)
volcano
## Warning: ggrepel: 61 unlabeled data points (too many overlaps). Consider increasing max.overlaps

Group - Intersection
intersection_results <- group_intersection(group, list_intersections, list_venns,
linda.output, fit_data,
raw_linda_results,
segment = segment,
level=level)
list_intersections <- intersection_results[[1]]
list_venns <- intersection_results[[2]]
venn <- intersection_results[[3]]
# show the results
venn

Basic statistics
uni_df <- merge(basic_univariate_statistics(linda_data,group),
raw_linda_results[[segment]][[group1]],
by="SeqID",all=TRUE)
## [1] "healthy"
## [1] "post_ltx"
## [1] "healthy"
## [1] "post_ltx"
uni_df[["final_sig"]] <- uni_df$SeqID %in% list_intersections[[paste(segment,level,comparison_name)]][["SeqID"]]
uni_statistics[[segment]][[paste(level,comparison_name)]] <- uni_df
# for comparison
new_name <- paste(level,comparison_name)
addWorksheet(supplements_wb, sheetName = new_name)
writeData(supplements_wb, sheet = new_name, uni_df, rowNames=FALSE)
Visualization
Heatmap visualizing the linDA’s logFoldChange for taxa with p <
0.1.
list_heatmap <- list_intersections[grep(paste(segment,level),
names(list_intersections),value=TRUE)]
p_heatmap_linda <- heatmap_linda(list_heatmap,colon_taxa_tab)
## Warning: The melt generic in data.table has been passed a data.frame and will attempt to redirect to the
## relevant reshape2 method; please note that reshape2 is superseded and is no longer actively developed, and
## this redirection is now deprecated. To continue using melt methods from reshape2 while both libraries are
## attached, e.g. melt.list, you can prepend the namespace, i.e. reshape2::melt(plot_df). In the next version,
## this warning will become an error.
## Warning: The melt generic in data.table has been passed a data.frame and will attempt to redirect to the
## relevant reshape2 method; please note that reshape2 is superseded and is no longer actively developed, and
## this redirection is now deprecated. To continue using melt methods from reshape2 while both libraries are
## attached, e.g. melt.list, you can prepend the namespace, i.e. reshape2::melt(plot_df). In the next version,
## this warning will become an error.
## Warning: The dcast generic in data.table has been passed a data.frame and will attempt to redirect to the
## relevant reshape2 method; please note that reshape2 is superseded and is no longer actively developed, and
## this redirection is now deprecated. Please do this redirection yourself like
## reshape2::dcast(plot_df_combined). In the next version, this warning will become an error.
p_heatmap_linda

Dot heatmap
dotheatmap_linda <- dot_heatmap_linda(list_heatmap,uni_df,colon_taxa_tab)
## Warning: The melt generic in data.table has been passed a data.frame and will attempt to redirect to the
## relevant reshape2 method; please note that reshape2 is superseded and is no longer actively developed, and
## this redirection is now deprecated. To continue using melt methods from reshape2 while both libraries are
## attached, e.g. melt.list, you can prepend the namespace, i.e. reshape2::melt(plot_df). In the next version,
## this warning will become an error.
## Warning: The melt generic in data.table has been passed a data.frame and will attempt to redirect to the
## relevant reshape2 method; please note that reshape2 is superseded and is no longer actively developed, and
## this redirection is now deprecated. To continue using melt methods from reshape2 while both libraries are
## attached, e.g. melt.list, you can prepend the namespace, i.e. reshape2::melt(plot_df). In the next version,
## this warning will become an error.
## Warning: The melt generic in data.table has been passed a data.frame and will attempt to redirect to the
## relevant reshape2 method; please note that reshape2 is superseded and is no longer actively developed, and
## this redirection is now deprecated. To continue using melt methods from reshape2 while both libraries are
## attached, e.g. melt.list, you can prepend the namespace, i.e. reshape2::melt(plot_df). In the next version,
## this warning will become an error.
## min_clr -0.9486582
## max_clr 6.236875
## min_log -4.029266
## max_log 3.847596
dotheatmap_linda

PSC effect
pre_LTx vs Healthy and Post_LTx vs Healthy
intersection
imagename
A <- list_intersections[[paste(segment,level,"healthy vs pre_ltx")]]$SeqID
B <- list_intersections[[paste(segment,level,"healthy vs post_ltx")]]$SeqID
C <- intersect(A,B)
psc_effect[[paste(segment,level)]] <- list_intersections[[paste(segment,level, "healthy vs pre_ltx")]][list_intersections[[paste(segment,level, "healthy vs pre_ltx")]]$SeqID %in% C,]
# see the results
psc_effect[[paste(segment,level)]]
Phylum level
level="phylum"
path_maaslin="../intermediate_files/maaslin/Q1_norway/Phylum/"
raw_linda_results_phylum[[segment]] <- list()
linda_results_phylum[[segment]] <- list()
Aggregate taxa
phylum_data <- aggregate_taxa(colon_asv_tab,
colon_taxa_tab,
taxonomic_level = "Phylum")
colon_phylum_tab <- phylum_data[[1]]
colon_phylum_taxa_tab <- phylum_data[[2]]
pre_ltx vs healthy
group <- c("healthy","pre_ltx")
comparison_name <- paste0(group[1], " vs ",group[2])
linDA
# prepare the data
linda_data <- binomial_prep(colon_phylum_tab,
colon_phylum_taxa_tab,
colon_metadata,
group, usage="linDA")
## Removing 1 ASV(s)
filt_colon_uni_data <- linda_data[[1]]
filt_colon_uni_taxa <- linda_data[[2]]
filt_colon_uni_metadata <- linda_data[[3]]
# fit the model
linda.obj <- linda(filt_colon_uni_data,
filt_colon_uni_metadata,
formula = '~ Group + (1|Patient)')
## 0 features are filtered!
## The filtered data has 274 samples and 8 features will be tested!
## Imputation approach is used.
## Fit linear mixed effects models ...
## Completed.
linda.output <- linda.obj$output
linda.output <- linda_renaming(linda.output, group)
group1 <- paste0(group[1], " vs ","Group",group[2])
for (grp in c(group1)){
raw_linda_results[[segment]][[grp]] <-
rawlinda.df(linda.output,
grp,
filt_colon_uni_data,
filt_colon_uni_taxa)
linda_results[[segment]][[grp]] <-
linda.df(linda.output,
grp,
filt_colon_uni_data,
filt_colon_uni_taxa)
}
# volcano plot
volcano_1 <- volcano_plot_linda(linda.output, group1,
taxa_table = filt_colon_uni_taxa) +
ggtitle(paste(group,collapse=" vs "))
## Using Phylum for naming
volcano <- ggarrange(volcano_1)
# see the plot
volcano

MaAsLin2
Volcano plot
volcano1 <- volcano_plot_maaslin(fit_data,filt_colon_uni_taxa) +
ggtitle(paste(group[1], "vs", group[2]))
volcano <- ggarrange(volcano1)
# see the results
volcano

Group - Intersection
intersection_results <- group_intersection(group, list_intersections, list_venns,
linda.output, fit_data,
raw_linda_results,
segment = segment,
level=level)
list_intersections <- intersection_results[[1]]
list_venns <- intersection_results[[2]]
venn <- intersection_results[[3]]
# show the results
venn

Basic statistics
uni_df <- merge(basic_univariate_statistics(linda_data,group),
raw_linda_results[[segment]][[group1]],
by="SeqID",all=TRUE)
## [1] "healthy"
## [1] "pre_ltx"
## [1] "healthy"
## [1] "pre_ltx"
uni_df[["final_sig"]] <- uni_df$SeqID %in% list_intersections[[paste(segment,level,comparison_name)]][["SeqID"]]
uni_statistics[[segment]][[paste(level,comparison_name)]] <- uni_df
# for comparison
new_name <- paste(level,comparison_name)
addWorksheet(supplements_wb, sheetName = new_name)
writeData(supplements_wb, sheet = new_name, uni_df, rowNames=FALSE)
pre_ltx vs post_ltx
group <- c("pre_ltx","post_ltx")
comparison_name <- paste0(group[1], " vs ",group[2])
linDA
# prepare the data
linda_data <- binomial_prep(colon_phylum_tab,
colon_phylum_taxa_tab,
colon_metadata,
group, usage="linDA")
## Removing 1 ASV(s)
## Removing 1 ASV(s)
filt_colon_uni_data <- linda_data[[1]]
filt_colon_uni_taxa <- linda_data[[2]]
filt_colon_uni_metadata <- linda_data[[3]]
# fit the model
linda.obj <- linda(filt_colon_uni_data,
filt_colon_uni_metadata,
formula = '~ Group + (1|Patient)')
## 0 features are filtered!
## The filtered data has 349 samples and 8 features will be tested!
## Pseudo-count approach is used.
## Fit linear mixed effects models ...
## Completed.
linda.output <- linda.obj$output
linda.output <- linda_renaming(linda.output, group)
group1 <- paste0(group[1], " vs ","Group",group[2])
for (grp in c(group1)){
raw_linda_results[[segment]][[grp]] <-
rawlinda.df(linda.output,
grp,
filt_colon_uni_data,
filt_colon_uni_taxa)
linda_results[[segment]][[grp]] <-
linda.df(linda.output,
grp,
filt_colon_uni_data,
filt_colon_uni_taxa)
}
# volcano plot
volcano_1 <- volcano_plot_linda(linda.output, group1,
taxa_table = filt_colon_uni_taxa) +
ggtitle(paste(group,collapse=" vs "))
## Using Phylum for naming
volcano <- ggarrange(volcano_1)
# see the plot
volcano

MaAsLin2
volcano1 <- volcano_plot_maaslin(fit_data,filt_colon_uni_taxa) +
ggtitle(paste(group[1], "vs", group[2]))
volcano <- ggarrange(volcano1)
# see the results
volcano

Group - Intersection
intersection_results <- group_intersection(group, list_intersections, list_venns,
linda.output, fit_data,
raw_linda_results,
segment=segment,
level=level)
list_intersections <- intersection_results[[1]]
list_venns <- intersection_results[[2]]
venn <- intersection_results[[3]]
# show the results
venn

Basic statistics
uni_df <- merge(basic_univariate_statistics(linda_data,group),
raw_linda_results[[segment]][[group1]],
by="SeqID",all=TRUE)
## [1] "pre_ltx"
## [1] "post_ltx"
## [1] "pre_ltx"
## [1] "post_ltx"
uni_df[["final_sig"]] <- uni_df$SeqID %in% list_intersections[[paste(segment,level,comparison_name)]][["SeqID"]]
uni_statistics[[segment]][[paste(level,comparison_name)]] <- uni_df
# for comparison
new_name <- paste(level,comparison_name)
addWorksheet(supplements_wb, sheetName = new_name)
writeData(supplements_wb, sheet = new_name, uni_df, rowNames=FALSE)
post_ltx vs healthy
group <- c("healthy","post_ltx")
comparison_name <- paste0(group[1], " vs ",group[2])
linDA
# prepare the data
linda_data <- binomial_prep(colon_phylum_tab,
colon_phylum_taxa_tab,
colon_metadata,
group, usage="linDA")
## Removing 6 ASV(s)
filt_colon_uni_data <- linda_data[[1]]
filt_colon_uni_taxa <- linda_data[[2]]
filt_colon_uni_metadata <- linda_data[[3]]
# fit the model
linda.obj <- linda(filt_colon_uni_data,
filt_colon_uni_metadata,
formula = '~ Group + (1|Patient)')
## 0 features are filtered!
## The filtered data has 207 samples and 8 features will be tested!
## Imputation approach is used.
## Fit linear mixed effects models ...
## Completed.
linda.output <- linda.obj$output
linda.output <- linda_renaming(linda.output, group)
group1 <- paste0(group[1], " vs ","Group",group[2])
for (grp in c(group1)){
raw_linda_results[[segment]][[grp]] <-
rawlinda.df(linda.output,
grp,
filt_colon_uni_data,
filt_colon_uni_taxa)
linda_results[[segment]][[grp]] <-
linda.df(linda.output,
grp,
filt_colon_uni_data,
filt_colon_uni_taxa)
}
# volcano plot
volcano_1 <- volcano_plot_linda(linda.output, group1,
taxa_table = filt_colon_uni_taxa) +
ggtitle(paste(group,collapse=" vs "))
## Using Phylum for naming
volcano <- ggarrange(volcano_1)
# see the plot
volcano

MaAsLin2
volcano1 <- volcano_plot_maaslin(fit_data,filt_colon_uni_taxa) +
ggtitle(paste(group[1], "vs", group[2]))
volcano <- ggarrange(volcano1)
volcano

Group - Intersection
intersection_results <- group_intersection(group, list_intersections, list_venns,
linda.output, fit_data,
raw_linda_results,
segment = segment,
level=level)
list_intersections <- intersection_results[[1]]
list_venns <- intersection_results[[2]]
venn <- intersection_results[[3]]
# show the results
venn

Basic statistics
uni_df <- merge(basic_univariate_statistics(linda_data,group),
raw_linda_results[[segment]][[group1]],
by="SeqID",all=TRUE)
## [1] "healthy"
## [1] "post_ltx"
## [1] "healthy"
## [1] "post_ltx"
uni_df[["final_sig"]] <- uni_df$SeqID %in% list_intersections[[paste(segment,level,comparison_name)]][["SeqID"]]
uni_statistics[[segment]][[paste(level,comparison_name)]] <- uni_df
# for comparison
new_name <- paste(level,comparison_name)
addWorksheet(supplements_wb, sheetName = new_name)
writeData(supplements_wb, sheet = new_name, uni_df, rowNames=FALSE)
Visualization
Heatmap visualizing the linDA’s logFoldChange for taxa with p <
0.1.
list_heatmap <- list_intersections[grep(paste(segment,level),
names(list_intersections),value=TRUE)]
p_heatmap_linda <- heatmap_linda(list_heatmap,colon_taxa_tab)
## Warning: The melt generic in data.table has been passed a data.frame and will attempt to redirect to the
## relevant reshape2 method; please note that reshape2 is superseded and is no longer actively developed, and
## this redirection is now deprecated. To continue using melt methods from reshape2 while both libraries are
## attached, e.g. melt.list, you can prepend the namespace, i.e. reshape2::melt(plot_df). In the next version,
## this warning will become an error.
## Warning: The melt generic in data.table has been passed a data.frame and will attempt to redirect to the
## relevant reshape2 method; please note that reshape2 is superseded and is no longer actively developed, and
## this redirection is now deprecated. To continue using melt methods from reshape2 while both libraries are
## attached, e.g. melt.list, you can prepend the namespace, i.e. reshape2::melt(plot_df). In the next version,
## this warning will become an error.
## Warning: The dcast generic in data.table has been passed a data.frame and will attempt to redirect to the
## relevant reshape2 method; please note that reshape2 is superseded and is no longer actively developed, and
## this redirection is now deprecated. Please do this redirection yourself like
## reshape2::dcast(plot_df_combined). In the next version, this warning will become an error.
p_heatmap_linda

Dot heatmap
dotheatmap_linda <- dot_heatmap_linda(list_heatmap,uni_df,colon_taxa_tab)
## Warning: The melt generic in data.table has been passed a data.frame and will attempt to redirect to the
## relevant reshape2 method; please note that reshape2 is superseded and is no longer actively developed, and
## this redirection is now deprecated. To continue using melt methods from reshape2 while both libraries are
## attached, e.g. melt.list, you can prepend the namespace, i.e. reshape2::melt(plot_df). In the next version,
## this warning will become an error.
## Warning: The melt generic in data.table has been passed a data.frame and will attempt to redirect to the
## relevant reshape2 method; please note that reshape2 is superseded and is no longer actively developed, and
## this redirection is now deprecated. To continue using melt methods from reshape2 while both libraries are
## attached, e.g. melt.list, you can prepend the namespace, i.e. reshape2::melt(plot_df). In the next version,
## this warning will become an error.
## Warning: The melt generic in data.table has been passed a data.frame and will attempt to redirect to the
## relevant reshape2 method; please note that reshape2 is superseded and is no longer actively developed, and
## this redirection is now deprecated. To continue using melt methods from reshape2 while both libraries are
## attached, e.g. melt.list, you can prepend the namespace, i.e. reshape2::melt(plot_df). In the next version,
## this warning will become an error.
## min_clr -4.405553
## max_clr 2.814299
## min_log -3.206076
## max_log 1.761042
dotheatmap_linda

PSC effect
pre_LTx vs Healthy and Post_LTx vs Healthy
intersection
imagename
A <- list_intersections[[paste(segment,level,"healthy vs pre_ltx")]]$SeqID
B <- list_intersections[[paste(segment,level,"healthy vs post_ltx")]]$SeqID
C <- intersect(A,B)
psc_effect[[paste(segment,level)]] <- list_intersections[[paste(segment,level, "healthy vs pre_ltx")]][list_intersections[[paste(segment,level, "healthy vs pre_ltx")]]$SeqID %in% C,]
# see the results
psc_effect[[paste(segment,level)]]